function a2i_out_sigma, ensemble

  nIt=n_elements(ensemble[0, *])
  dOut=fltarr(2, n_elements(ensemble[*, 0]))
  for ti=0, n_elements(ensemble[*, 0])-1 do begin
    sorted=reform(ensemble[ti, *])
    sorted=sorted[sort(sorted)]
    dOut[0, ti]=sorted[0.16*nIt]
    dOut[1, ti]=sorted[0.84*nIt]
  endfor
  
  return, dOut
  
end

function a2i_out_kz, inTime, in, inEnsemble, average=average, growth=growth

  window=0.7
  kzIt=4

  if keyword_set(average) then begin
    averageScaling=4.
  endif else begin
    averageScaling=1.
  endelse

  nIt=n_elements(inEnsemble[0, 0, *])
  
  if size(reform(in), /n_dimensions) eq 2 then begin
    toSmooth=total(in, 1)/averageScaling
  endif else begin
    toSmooth=in
  endelse

  kzOut=mr_kz_filter(inTime, toSmooth, window=window, Iterations=kzIt, growth=growth)
  outTime=kzOut.x
  out=kzOut.Y

  outEnsemble=fltarr(n_elements(inTime) - keyword_set(growth), nIt)  
  for it=0, nIt-1 do begin
    if size(reform(in), /n_dimensions) eq 2 then begin
      toSmooth=total(reform(inEnsemble[*, *, it]), 1)/averageScaling
    endif else begin
      toSmooth=inEnsemble[*, it]
    endelse
    kzOut=mr_kz_filter(inTime, toSmooth, window=window, Iterations=kzIt, growth=growth)
    outEnsemble[*, it]=kzOut.Y
  endfor

  ;Calculate uncertainty
  dOut=a2i_out_sigma(outEnsemble)

  return, {time:outTime, average:out, sigma:dOut}

end

function a2i_out_annual_function, in
  
  case size(in, /n_dimensions) of
    1: out=fltarr(n_elements(in)/12)
    2: out=fltarr(n_elements(in[0, *])/12)
  endcase
  for ti=0, n_elements(out)-1 do begin
    case size(in, /n_dimensions) of
      1: out[ti]=mean(in[ti*12:(ti+1)*12-1], /nan)
      2: out[ti]=total(in[*, ti*12:(ti+1)*12-1], /nan)
    endcase
  endfor
  return, out
  
end

function a2i_out_annual, inTime, in, inEnsemble, average=average

  nYears=n_elements(inTime)/12
  nIt=n_elements(inEnsemble[0, 0, *])
  
  outTime=fltarr(nYears)
  out=fltarr(nYears)
  outSigma=fltarr(nYears)
  outEnsemble=fltarr(nYears, nIt)

  outTime=a2i_out_annual_function(inTime)
  out=a2i_out_annual_function(in)
  for it=0, nIt-1 do begin
    outEnsemble[*, it]=a2i_out_annual_function(reform(inEnsemble[*, *, it]))
  endfor

  if keyword_set(average) then begin
    out=out/12./4.
    outEnsemble=outEnsemble/12./4.
  endif else begin
    out=out/12.
    outEnsemble=outEnsemble/12.
  endelse

  ;Calculate uncertainty
  dOut=a2i_out_sigma(outEnsemble)

  return, {time:outTime, average:out, sigma:dOut}

end

function a2i_out_monthly, inTime, in, inEnsemble, average=average

  timesize=n_elements(inTime)
  nIt=n_elements(inEnsemble[0, 0, *])
  
  outTime=fltarr(timeSize)
  out=fltarr(timeSize)
  outSigma=fltarr(timeSize)
  outEnsemble=fltarr(timeSize, nIt)
  
  outTime=inTime
  out=total(in, 1)
  for it=0, nIt-1 do begin
    outEnsemble[*, it]=total(reform(inEnsemble[*, *, it]), 1)
  endfor
  
  if keyword_set(average) then begin
    out=out/4.
    outEnsemble=outEnsemble/4.
  endif
  
  ;Calculate uncertainty
  dOut=a2i_out_sigma(outEnsemble)

  return, {time:outTime, average:out, sigma:dOut}
  
end


pro a2i_out_write_semi_hemispheric, case_name, folder, pollutant, description, data, $
  prefix=prefix_in, pollutant_string=pollutant_string, unit=unit

  if keyword_set(prefix_in) then begin
    prefix=prefix_in
  endif else begin
    prefix=''
  endelse

  file_str=a2i_filestr(/input,  + case_name + '/output/' + folder + '/' + prefix + pollutant)

  ;Save IDL variables
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  time=data.time
  average=data.average
  sigma=data.sigma
  
  save, filename=file_str + '.sav', time, average, sigma, pollutant, pollutant_string, unit

  ;Write CSV files
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

  openw, fun, file_str + '.csv', /get_lun

    if folder eq 'emissions' then begin
      printf, fun, '%' + pollutant + ' ' + description + ',,,,,,'
      printf, fun, '%' + 'Matt Rigby University of Bristol (matt.rigby@bristol.ac.uk),,,,,,'
      printf, fun, '%' + systime() + ',,,,,,'
      printf, fun, 'Time, Global, 30N-90N, 0N-30N, 30S-0S, 90S-30S,'
      for ti=0, n_elements(data.time)-1 do begin
        printf, fun, time[ti], total(average[*, ti]), average[*, ti], format="(6(f16.5, ','))"
      endfor
    endif else begin
      printf, fun, '%' + pollutant + ' ' + description + ',,,,,'
      printf, fun, '%' + 'Matt Rigby University of Bristol (matt.rigby@bristol.ac.uk),,,,,'
      printf, fun, '%' + systime() + ',,,,,'
      printf, fun, 'Time, 30N-90N, 0N-30N, 30S-0S, 90S-30S,'
      for ti=0, n_elements(data.time)-1 do begin
        printf, fun, time[ti], average[*, ti], format="(5(f16.5, ','))"
      endfor
    endelse
    
  close, fun
  free_lun, fun
  
end

pro a2i_out_write_global, case_name, folder, pollutant, description, data, $
  prefix=prefix_in, pollutant_string=pollutant_string, unit=unit, column_header=column_header

  if keyword_set(prefix_in) then begin
    prefix=prefix_in
  endif else begin
    prefix=''
  endelse

  if ~keyword_set(column_header) then begin
    column_header=description
  endif
  
  file_str=a2i_filestr(/input,  + case_name + '/output/' + folder + '/' + prefix + pollutant)
  
  ;Save IDL variables
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  time=data.time
  average=data.average
  sigma=data.sigma
  
  save, filename=file_str + '.sav', time, average, sigma, pollutant, pollutant_string, unit
  
  ;Write CSV files
  ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
  
  openw, fun, file_str + '.csv', /get_lun
  
    printf, fun, '%' + pollutant + ' ' + description + ',,,,'
    printf, fun, '%' + 'Matt Rigby University of Bristol (matt.rigby@bristol.ac.uk),,,,'
    printf, fun, '%' + systime() + ',,,,'
    printf, fun, 'Time, ' + column_header + ', 16%ile, 84%ile,'
    
    for ti=0, n_elements(data.time)-1 do begin
      printf, fun, time[ti], average[ti], sigma[0, ti], sigma[1, ti], format="(4(f16.5, ','))"
    endfor
  
  close, fun
  free_lun, fun
  
end


pro a2i_out, case_name, nIt=nIt

  if n_elements(case_name) eq 0 then begin
    print, 'Case_name needed'
    return
  endif

  seed=systime(/seconds)

  restore, a2i_filestr(/input, case_name + '/output/output.sav'), /relaxed

  if n_elements(H) eq 0 then begin
    H=fltarr(x['NSTATE'], y['NMEASUREMENTS'])
    P=fltarr(x['NSTATE'], x['NSTATE'])
    H[temporary(Hwh)]=temporary(Hsp)
    P[temporary(Pwh)]=temporary(Psp)
  endif

  if ~keyword_set(nIt) then begin
    nIt=100
  endif

  ;Get lifetimes
  lifetime=!null
  openr, fun, a2i_filestr(/input, case_name + '/output/lifetimes.csv'), /get_lun
  str=''
  readf, fun, str
  for pi=0, parameters['NPOLLUTANTS']-1 do begin
    str=''
    readf, fun, str
    data=strSplit(str, ',', /extract)
    lifetime=[lifetime, float(data[2])]
  endfor
  close, fun
  free_lun, fun


  for pi=0, parameters['NPOLLUTANTS']-1 do begin
    
    whY=where(y['POLLUTANT'] eq pi, cYS)
    whX=where(x['POLLUTANT'] eq pi or x['POLLUTANT'] lt 0, cXS)
    Hs=H[*, whY]
    Hs=Hs[whX, *]
    Ps=P[*, whX]
    Ps=Ps[whX, *]
    xS=x['X', whX]
    xIDS=x['ID', whX]
    xBoxS=x['BOX', whX]
    yMS=y['MODEL', whY]
    yS=y['Y', whY]
    dyS=y['ERROR', whY]
    yBoxS=y['BOX', whY]
    yTiS=y['TI', whY]

    timeS=parameters['STARTY'] + yTiS[where(yBoxS eq 0)]/12. + 1./24.

    ;Averages
    ; mole fractions
    mfS=fltarr(4, cYS/4)
    dmfS=fltarr(4, cYS/4)
    mfMS=fltarr(4, cYS/4)
    for bi=0, 3 do begin
      wh=where(yBoxS eq bi)
      mfS[bi, *]=yS[wh]
      dmfS[bi, *]=dyS[wh]
      mfMS[bi, *]=yMS[wh]
    endfor

    ; emissions
    qS=model['EMISSIONS', pi, *, *]

    ;Calculate ensembles
    if total(abs(PS)) gt 0. then begin
      xSEn=mr_cholesky_covariance(xS, PS, nIt)
    endif else begin
      xSEn=fltarr(cXS, nIt)
      for it=0, nIt-1 do begin
        xSEn[*, it]=xS
      endfor
    endelse

    ySEn=fltarr(cYS, nIt)
    mfMSEn=fltarr(4, cYS/4, nIt)
    qSEn=fltarr(n_elements(model['EMISSIONS', pi, *, 0]), n_elements(model['EMISSIONS', pi, 0, *]), nIt)
    
    ;Calculate burden
    mfAn=mean(mfmS, dimension=1)
;    mfAn=rebin(mfAn, n_elements(mfAn)/12)
    case parameters['UNIT', pi] of
      'ppq': begin
        scale_factor=1.e-15
        emissions_scale_factor=1.e6
      end
      'ppt': begin
        scale_factor=1.e-12
        emissions_scale_factor=1.e9
      end
      'ppb': begin
        scale_factor=1.e-9
        emissions_scale_factor=1.e12
      end
      'ppm': begin
        scale_factor=1.e-6
        emissions_scale_factor=1.e15
      end
    endcase
    burdenS=mfAn*scale_factor*!mr_const.mAtm/!mr_const.mAir*1000.*parameters['MOL_MASS', pi] ;g

    if lifetime[pi] gt 0. and lifetime[pi] lt 1.e4 then begin
      lifetime_error_fractional=$
        (parameters['GLOBAL_LIFETIME_ERROR', pi]/lifetime[pi]*(burdenS/emissions_scale_factor))/$
        total(qS, 1)
    endif else begin
      ;MINIMUM 1% uncertainty (model transport error)
      lifetime_error_fractional=replicate(0.01, n_elements(qS[0, *]))
    endelse


    ;Generate ensembles of mole fractions and emissions
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    for n=0, nIt-1 do begin
      
      ;Mole fraction
      ySEn[*, n]=yMS + Hs##(xSEn[*, n] - xS)
      for bi=0, 3 do begin
        wh=where(yBoxS eq bi)
        mfMSEn[bi, *, n]=ySEn[wh, n]
      endfor
      ; scale error
      mfMSEn[*, *, n]=mfMSEn[*, *, n]*(1. + parameters['SCALE_ERROR_PLUS', pi]*randomn(seed))[0]

      ;Emissions
;      qSEn[*, *, n]=model['EMISSIONS', pi, *, *]
      for bi=0, 3 do begin
        wh=where(xBoxS eq bi and xIDS eq 'Q', count)
        xTemp=reform(xSEn[wh, n] - xS[wh])
        xTemp=congrid(xTemp, n_elements(model['EMISSIONS', pi, 0, *]), /center)
        qSEn[bi, *, n]=model['EMISSIONS', pi, bi, *]*(1. + xTemp)
      endfor

      ; scale error and lifetime error
      qSEN[*, *, n]=qSEN[*, *, n]*(1. + $
        replicate(parameters['SCALE_ERROR_PLUS', pi]*randomn(seed), 4, n_elements(qSEN[0, *, 0])) + $
        transpose(rebin(lifetime_error_fractional*randomn(seed), n_elements(qSEN[0, *, 0]), 4)))

    endfor

    ;Calculate output quantities
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    outMFmonth=a2i_out_monthly(timeS, mfMS, mfMSEn, /average)

    outMFan=a2i_out_annual(timeS, mfMS, mfMSEn, /average)
    outQan=a2i_out_annual(timeS, qS, qSEn)

    qGmonth=fltarr((size(qS))[1:(size(qS))[0]])
    qGmonthEN=fltarr((size(qSEN))[1:(size(qSEN))[0]])
    qGmonth[*, 12:-1]=qS[*, 12:-1] - qS[*, 0:-12]
    for n=0, nIt-1 do qGmonthEN[*, 12:-1, n]=qSEN[*, 12:-1, n]-qSEN[*, 0:-12, n]
    outQGan=a2i_out_annual(timeS, qGmonth, qGmonthEN)

    nrebin=n_elements(timeS)/12
    qGRmonth=fltarr(nrebin)
    qGRmonthEN=fltarr(nrebin, nIt)
    qGRmonth[1:-1]=rebin(total(qS[*, 12:-1], 1), nrebin-1)/rebin(total(qS[*, 0:-13], 1), nrebin-1)-1.
    for n=0, nIt-1 do qGRmonthEN[1:-1, n]=rebin(total(qSEN[*, 12:-1, n], 1), nRebin-1)/$
      rebin(total(qSEN[*, 0:-13, n], 1), nRebin-1)-1.
    dqGRmonth=a2i_out_sigma(qGRmonthEN)
    outQGRan={time:rebin(timeS-0.5, n_elements(qGRmonth)), average:qGRmonth, sigma:dqGRmonth}

;    ;Percent emissions change relative to 2000 (change this to be a year that you input)
;    dummy=min(abs(timeS - 1999.5), tiRef)
;    qPC=qS / rebin(mean(qS[*, tiRef + 0:tiRef + 11], dim=2), 4, n_elements(qS[0, *]))*100. - 100.
;    qPCEn=qSEn
;    for n=0, nIt-1 do begin
;      qPCEn[*, *, n] = qSEn[*, *, n] / rebin(mean(qSEn[*, tiRef + 0:tiRef + 11, n], dim=2), 4, n_elements(qSEn[0, *, 0]))*100.-100.
;    endfor
;    outQPC=a2i_out_kz(timeS, qPC, qPCEn, /average)


    outMFKZ=a2i_out_kz(timeS, mfMS, mfMSEn, /average)
    outQKZ=a2i_out_kz(timeS, qS, qSEn)

    outMFGKZ=a2i_out_kz(timeS, mfMS, mfMSEn, /average, /growth)
    outQGKZ=a2i_out_kz(timeS, qS, qSEn, /growth)

    ;semi-hemispheric growth
    shaverage=!null
    shsigma=!null
    shtime=!null
    for bi=0, 3 do begin
      out=a2i_out_kz(timeS, reform(mfMS[bi, *]), reform(mfMSEn[bi, *, *]), /growth)
      shaverage=[[shaverage], [out.average]]
      shsigma=[[shsigma], [out.sigma]]
    endfor
    outMFGSHKZ={average:transpose(shaverage), sigma:transpose(shsigma), time:out.time}


    ;Write outputs
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
    
    case parameters['UNIT', pi] of
      'ppq': emissions_unit='Mg/yr'
      'ppt': emissions_unit='Gg/yr'
      'ppb': emissions_unit='Tg/yr'
      'ppm': emissions_unit='Pg/yr'
    endcase

    ;Semi-hemispheric modeled mole fractions
    a2i_out_write_semi_hemispheric, case_name, $
      'mf', $
      parameters['POLLUTANTS', pi], $
      'Modeled mole fraction (' + parameters['UNIT', pi] + ') output from AGAGE 12-box model', $
      {time:timeS, average:mfMS, sigma:fltarr(4, n_elements(timeS))}, $
      prefix='mf_model_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=parameters['UNIT', pi]

    ;Semi-hemispheric mole fractions
    a2i_out_write_semi_hemispheric, case_name, $
      'mf', $
      parameters['POLLUTANTS', pi], $
      'AGAGE mole fraction (' + parameters['UNIT', pi] + ') averaged into semi-hemispheres', $
      {time:timeS, average:mfS, sigma:transpose([[[mfS - dmfS]], [[mfS + dmFS]]], [2, 0, 1])}, $
      prefix='mf_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=parameters['UNIT', pi]

    ;Global annual mole fraction
    a2i_out_write_global, case_name, $
      'mf', $
      parameters['POLLUTANTS', pi], $
      'Global mole fraction (' + parameters['UNIT', pi] + ') output from AGAGE 12-box model', $
      outMFan, $
      prefix='mf_global_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=parameters['UNIT', pi]

    ;Global monthly mole fraction
    a2i_out_write_global, case_name, $
      'mf', $
      parameters['POLLUTANTS', pi], $
      'Global mole fraction (' + parameters['UNIT', pi] + ') output from AGAGE 12-box model', $
      outMFmonth, $
      prefix='mf_global_monthly_', pollutant_string=parameters['POLLUTANT_STRING', pi], $
      unit=parameters['UNIT', pi]

    ;Global smoothed growth rate
    a2i_out_write_global, case_name, $
      'mf', $
      parameters['POLLUTANTS', pi], $
      'Global growth rate (' + parameters['UNIT', pi] + '/yr) output from AGAGE 12-box model', $
      outMFGKZ, $
      prefix='mf_growth_', pollutant_string=parameters['POLLUTANT_STRING', pi], $
      unit=parameters['UNIT', pi] + '/yr'

    ;Semi-hemispheric growth rate
    a2i_out_write_semi_hemispheric, case_name, $
      'mf', $
      parameters['POLLUTANTS', pi], $
      'Modeled mole fraction growth rate (' + parameters['UNIT', pi] + '/yr) output from AGAGE 12-box model', $
      outMFGSHKZ, $
      prefix='mf_growth_sh_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=parameters['UNIT', pi] + '/yr'

    ;Annual radiative forcing
    a2i_out_write_global, case_name, $
      'mf', $
      parameters['POLLUTANTS', pi], $
      'Radiative forcing (mW/m2) output from AGAGE 12-box model', $
      {time:outMFAn.time, $
        average:mr_radiative_forcing(outMFAn.average, parameters['POLLUTANTS', pi]), $
        sigma:mr_radiative_forcing(outMFAn.sigma, parameters['POLLUTANTS', pi])}, $
      prefix='rf_global_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=parameters['UNIT', pi]


    ;EMISSIONS
    ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

    ;Global annual emissions
    a2i_out_write_global, case_name, $
      'emissions', $
      parameters['POLLUTANTS', pi], $
      'Global emissions (' + emissions_unit + ') output from AGAGE 12-box model. Assumed lifetime = ' + $
      string(lifetime[pi], format='(f10.1)') + ' years', $
      outQan, $
      prefix='emissions_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=emissions_unit, $
      column_header='Global emissions (' + emissions_unit + ')'

    ;Global annual emissions growth rate
    a2i_out_write_global, case_name, $
      'emissions', $
      parameters['POLLUTANTS', pi], $
      'Global emissions (' + emissions_unit + '/yr) output from AGAGE 12-box model. Assumed lifetime = ' + $
        string(lifetime[pi], format='(f10.1)') + ' years', $
      outQGan, $
      prefix='emissions_an_growth_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=emissions_unit + '/yr', $
      column_header='Global emissions (' + emissions_unit + '/yr)'

    ;Global annual emissions growth rate ratio
    a2i_out_write_global, case_name, $
      'emissions', $
      parameters['POLLUTANTS', pi], $
      'Global emissions (kg/kg/yr) output from AGAGE 12-box model. Assumed lifetime = ' + $
      string(lifetime[pi], format='(f10.1)') + ' years', $
      outQGRan, $
      prefix='emissions_an_growth_ratio_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit='kg/kg/yr', $
      column_header='Global emissions (kg/kg/yr)'
    
    ;Global emissions growth rate
    a2i_out_write_global, case_name, $
      'emissions', $
      parameters['POLLUTANTS', pi], $
      'Global emissions growth rate (' + emissions_unit + '/yr) output from AGAGE 12-box model' + $
        string(lifetime[pi], format='(f10.1)') + ' years', $
      outQGKZ, $
      prefix='emissions_growth_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=emissions_unit + '/yr', $
      column_header='Global emissions growth (' + emissions_unit + '/yr)'

    ;Semi-hemispheric emissions
    a2i_out_write_semi_hemispheric, case_name, $
      'emissions', $
      parameters['POLLUTANTS', pi], $
      'Emissions (' + emissions_unit + ') output from AGAGE 12-box model. Assumed lifetime = ' + $
      string(lifetime[pi], format='(f10.1)') + ' years', $
      {time:timeS, average:qS, sigma:fltarr(4, n_elements(timeS))}, $
      prefix='emissions_sh_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=emissions_unit

;    ;Global Percentage emissions
;    a2i_out_write_global, case_name, $
;      'emissions', $
;      parameters['POLLUTANTS', pi], $
;      'Global emissions percent change relative to 2000 (%) output from AGAGE 12-box model. Assumed lifetime = ' + $
;      string(lifetime[pi], format='(f10.1)') + ' years', $
;      outQPC, $
;      prefix='emissions_percent_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=emissions_unit + '/yr', $
;      column_header='Global emissions relative to 2000 (%)'

    ;Global annual CO2-equivalent emissions
    a2i_out_write_global, case_name, $
      'emissions', $
      parameters['POLLUTANTS', pi], $
      'Global CO2-equivalent emissions (' + emissions_unit + ') output from AGAGE 12-box model. Assumed lifetime = ' + $
      string(lifetime[pi], format='(f10.1)') + ' years', $
      {time:outQAn.time, $
        average:mr_co2e(outQAn.average, parameters['POLLUTANTS', pi]), $
        sigma:mr_co2e(outQAn.sigma, parameters['POLLUTANTS', pi])}, $
      prefix='emissions_co2e_', pollutant_string=parameters['POLLUTANT_STRING', pi], unit=emissions_unit

    print, '... done ' + parameters['POLLUTANTS', pi]

  endfor


end